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Abstract 

We formulate a time-dependent density functional theory (TDDFT) in terms of the density 
matrix to study ultrafast phenomena in semiconductor structures. A system of equations for the 
density matrix components, which is equivalent to the time-dependent Kohn-Sham equation, is 
derived. From this we obtain a TDDFT version of the semiconductor Bloch equations, where the 
electronic many-body effects are taken into account in principle exactly. As an example, we study 
the optical response of a three-dimensional two-band insulator to an external short-time pulsed 
laser field. We show that the optical absorption spectrum acquires excitonic features when the 
exchange-correlation potential contains a 1/q 2 Coulomb singularity. A qualitative comparison of 
the TDDFT optical absorption spectra with the corresponding results obtained within the Hartree- 
Fock approximation is made. 

PACS numbers: 71.10.-w, 71.15.Mb, 71.45. Gm 
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I. INTRODUCTION 



Due to the demands of modern electronics, semiconductor devices are becoming smaller 
and faster, which means that applied external fields cause strongly inhomogeneous and 
nonequilibrium processes in such systems. A powerful approach to study dynamical proper- 
ties of potentially useful devices is to apply short (femto- or picosecond) electric field pulses 
and to measure the response of the system. Recent progress in ultrafast laser pulse exper- 
imental techniques allows one to study the physical processes in such systems with a very 
high precision (for a review, see e.g. Refs. Il|j2j). For example, with this technique one can 
measure nonequilibrium energy and momentum distributions, or the dynamics of excited 
states. Therefore, it is necessary to develop theoretical tools to describe these experiments, 
as well as to understand ultrafast processes in small semiconductor devices in general. 

The theoretical description of ultrafast phenomena triggered by short-pulse laser fields is a 
complicated problem due to several reasons. One of the most difficult tasks is to take into ac- 
count many-particle correlation effects properly. An external pulse field causes the following 
main effects in the system: i) direct electron photoemission; ii) inverse electron photoemis- 
sion, and iii) absorption processes. The first two phenomena can be described in terms of free 
quasi-particles, which makes the problem relatively simple. A proper description of optical 
absorption spectra is a much more complicated task due to quasi-particle correlation effects. 
In particular, external laser pulses can create excitons, or coupled electron-hole pairs. The 
problem of correctly describing optical absorption spectra with excitonic features in the case 
of applied short-time laser pulses is one of the great challenges in condensed matter physics. 

For small devices in the presence of a short-pulse field typical time scales are shorter than 
the Coulomb scattering time (see, for example, Ref. O), which means that one cannot treat 
the Coulomb interaction effects by using a simple Boltzmann equation approach, where all 
the Coulomb effects are "hidden" in a scattering time parameter r. Similarly, the semicon- 
ductor Bloch equation (SBE) approach,- based on the Hartree-Fock (HF) approximation, 
and other mean-field methods have difficulties under these circumstances, because of the 
presence of strong fluctuations. In principle, Coulomb interaction effects can be taken into 
account in a systematic way by using non-equilibrium Green function techniques,- 1 ^ similar 
to the equilibrium case. Unfortunately, this technique becomes numerically very compli- 
cated in a strongly nonequilibrium situation, since in this case the Green functions depend 
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on two or more time arguments, which puts high requirements on the computer memory 
size and makes the numerical analysis very time-consuming.-^ 

In this paper, we discuss an alternative and potentially very powerful approach to study 
these kinds of problems, based on density-functional theory (DFT),- and in particular 
its time-dependent generalization (TDDFT).^^* 1 ! In this approach, numerical calculations 
should be much less time-consuming compared to the Green's function method. TDDFT 
has been successfully applied to describe molecular excitations;— however, this approach 
has some difficulties in describing extended systems.— It is known that the standard local- 
density approximation (LDA) and generalized gradient approximation (GGA) for the DFT 
exchange-correlation (xc) potential cannot be applied to describe some effects beyond the 
ground state properties in extended systems, like the energy band gaps and excitonic effects 
in the optical absorption spectra. Therefore, in order to apply TDDFT to study ultra- 
fast processes, and in particular to describe correctly the optical absorption spectra in such 
systems, it is necessary to find suitable xc potentials. 

For weak and smooth external fields, such a potential can constructed by using the 
many-body Bethe-Salpeter equation (BSE) approacbr^^. In fact, in this case TDDFT 
calculations with an xc potential extracted from the BSE give very good results for the 
optical absorption spectra in a set of bulk semiconductors. Unfortunately, the BSE-TDDFT 
approach cannot be used directly to construct an xc potential for excitation with strong short 
pulses, since here the linear response theory cannot be applied, and the equations in time 
domain depend on many time variables, which makes numerical solution extremely difficult 
to obtain (similar to non-equilibrium Green's functions). Therefore, it would be extremely 
useful to find a simple xc potential which will allow one to use TDDFT straightforwardly 



to study the optical response of a system in the time domain. It was shown in Refs. 



13, 



18l that for smooth and weak fields such a potential exists. Namely, the exact-exchange 



approximation for the time-dependent optimized effective potential (XX-TDOEP) results 
in optical absorption spectra with pronounced excitonic features, by solving the problem in 
the frequency domain. Kim and Gorling^ 1 ^ showed that the main reason for this was the 
presence of 1/q 2 Coulomb singularity in the exchange energy kernel. The XX-TDOEP has 
recently been extended into the nonlinear, real-time domain for simple quasi-one- dimensional 
quantum well systems.— 

The purpose of this paper is twofold: 1) we will formulate a general TDDFT approach 
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to study optical interband excitations in terms of the Kohn-Sham density matrix. The 
corresponding system of equations has the formal simplicity of the SBEs, which allows one 
to solve it directly in the time domain, contrary to the many-body Green function approach, 
where it is difficult to treat the problem numerically. At the same time, this approach has 
a great formal advantage in comparison with the SBE formalism, since in TDDFT the 
Coulomb interaction effects are treated in principle exactly. 2) We will show that by using 
simple exchange-only functionals, the essence of excitonic features can be captured in a 
relatively simple manner. As an example, we consider the model of a three-dimensional 
two-band bulk insulator for different local xc potentials and show that its optical absorption 
spectra contain qualitatively correct excitonic features, whenever the xc energy kernel has a 
1/q 2 singularity. 

The paper is organized as follows: We introduce a general TDDFT-SBE formalism in 
terms of the Kohn-Sham density matrix in Section [III In Section lITTj we derive the TDDFT- 
SBE formalism and in Section [IV] we apply this formalism to study the optical absorption 
spectra for a three-dimensional two-band model insulator with different xc potentials and 
compare the results with HF. Conclusions are presented in Section |V] Some technical details 
are given in Appendix. 

We use Hartree atomic units (e 2 = m = h = 1) throughout the paper. 



II. GENERAL FORMALISM 

The general DFT Hamiltonian for a many-electron system in a solid (in the Born- 
Oppenheimer approximation) can be written in the following form: 

H = -^- + V nud (r) + V H [n){v) + V xc [n)(r), (1) 

where V nuc i(r) is the nuclear potential for the electrons and 

V H [n]{v) = J dv'^-^ (2) 

is the Hartree potential, where n(r) is the density of electrons. All many-body effects 
beyond Hartree are described by the scalar xc potential V^ c [n](r). In particular, in the LDA 
exchange-only case: 

/o\l/3 

V xLDA (n(r)) = -(-) n^r). (3) 
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In order to describe the ground-state properties of the system governed by the Hamilto- 
nian Eq. (JT]), one solves the stationary Kohn-Sham (KS) equation 

tf(r)4*(r)=4^(r), (4) 

and find bands U with spectra e£, where k is the crystal momentum and l, L = Vi, Ci(i = 1,2,...) 
are the labels for the valence (v) and conduction (c) bands. The electron density can be 
found self-consistently: 

n(r)=2j2\^(r)\ 2 e(e F -eZ), (5) 

i,k 

where £p is the Fermi energy and the summation is performed over the occupied (valence) 
band states. Equation (jSJ) is valid in the case of zero temperature, which we consider in this 
paper. Finite temperatures would require introduction of the Fermi distribution function 
into Eq. (JSI). 

In order to study the nonequilibrium case when an external electric field E(i) is switched 
on at time t = to, the Hamiltonian Eq. (TjQ) must be modified in the standard way: an 
external electromagnetic vector potential A ext (r,t) should be added by making the usual 
substitution V — > V — (i/c)A ext (r,t) and by adding a scalar potential term <p ext (r,t) to the 
Hamiltonian. The electric field is connected with (p ex t(r,t) and A ext (r, t) in the following 
way: 

E(r, t) = -V^,(r, t) - - f^^ - (6) 

For an extended system one needs to preserve the periodicity, therefore the vector external 
potential must be used. This makes the problem technically more complicated in comparison 
to scalar potentials. However, it can be shown that in situations when the characteristic 
field frequency is bigger than the level spa cing , one can work with scalar potentials instead 



of vector potentials (see, for example, Ref. |2Q| ). Therefore, we shall consider the case when 



A ext (r, t) = 0, and for simplicity assume that the electric field is space- independent. Strictly 
speaking, for external homogeneous fields one needs to use current-TDDFT in order to study 
the response of the system.— 1 ^ One then gets the macroscopic current, which allows one to 
satisfy the periodicity condition. In this paper, however, we use the usual approximation 4 
where the external homogeneous electric field is described by the following scalar potential: 

y> est (r,t) = -E(t)r. (7) 
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Due to the presence of an electric field, the Hamiltonian is explicitly time-dependent. 
Moreover, the Hartree and xc potential terms in the Hamiltonian become also time- 
dependent. Such a time-dependent problem is described by the time-dependent KS equation: 

zJU(r,t) = #(r,t)*(r,t). (8) 

In the adiabatic approximation, which we use in this paper, the xc potential can be formally 
obtained from the stationary potentials by assuming that the electron density becomes 
time-dependent, n(r) — > n(r,t). Eq. (jSJ) must be solved self-consistently together with the 
corresponding time-dependent particle number equation (J5]). Since the initially occupied 
states are the valence states r/>£* (r), described by the band index v l and the momentum 
k, the evolution of the system can be completely described by the time-evolution of these 
states, i.e. by finding the corresponding time-dependent wave functions ^(r, i), such that 
(r, t ) = z/>k (r), from Eq. (jBJ) and a time-dependent generalization of Eq. (jSJ): 

n(r,t) = 2j2\n(r,t)\ 2 6(sF-eZ). (9) 

i,k 

In the following, we express the time-dependent wave functions as linear combinations of 
the ground-state wave functions: 

n (r, *) = £ [eg* (* W (r) + (t)^ (r)] , (10) 



where (i) are momentum and time-dependent complex coefficients, which satisfy the 
following initial condition: c k y(to) = ^kq^^e - *^ °. The time evolution of the system 
can thus be found by determining the coefficients '(£) and c^{t) for Eq. (TlTJT) . How- 
ever, in order to solve the problem, it is more convenient to introduce the density matrix: 
Pk-qp ™ (*) = c kq m (*) c kp" (*) > which is useful for defining physical quantities like occupa- 
tion of the states and optical transitions (see next Section). TDDFT in the density-matrix 
representation is a method which allows one to solve the problem in principle exactly, since 
it is an exact reformulation of the time- dependent Kohn-Sham formalism. This method was 
already introduced to study intersubband processes in quantum wells^ and the electrical 
conductivity in dissipative models of molecular devices^ (see also Ref. 



25 



where a similar 



approach was discussed). Here we develop a technique which can be applied in more general 
cases with a continuous electron spectrum, including interband transitions in solids. 
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The density matrix satisfies the following equation of motion: 

^p^Ht) = mt),p}^ = £ [H l ^(t) P ^ l ;(t) - P ^hm 1 ^ (t)} , (ii) 

where the Hamiltonian matrix elements Hj^ n (t) are 

HtjHt)= I dnfc*{r)H (i)<«(r) = efrS^ + EWd^ +V%£{t) + v2£(t), (12) 

J cell 

and the space integration is performed over a unit cell. In Eq. ffl2l) . 

= / rfr<™*(r)r^(r) (13) 

J cell 



I I I I 

are the dipole matrix elements, and V^™(£) and V^"(t) are the matrix elements for the 
difference between the time-dependent and ground state (at t < t ) Hartree and xc poten- 
tials: 

V%£{t)= I dr^t*(r)[V H [n](r,t)-V H [n}(v,t )}^( r ), (14) 

J cell 



I I' 

and similar for V^"(t). The density matrix satisfies the following initial condition: 

Pk-qpH*o) = S kcl 8 k p8 Vi i m 8 v . Vn , (15) 

which corresponds to the situation where all states in the valence bands are initially occupied. 
The particle density has the following form in terms of the density matrix elements: 

n(r,f) = 2 J2 Plif^W^^F - e%). (16) 

j,Z m ,/^,k,q,p 

Solution of the Liouville-von Neumann equation (fTTl) allows one to study physical properties 
of the system. In particular, one can find the polarization as 

D (*)= E Pki C Wd^. (17) 

Below, we shall use this formalism to study the optical response of a three-dimensional 
two-band insulator, by solving Eqs. fTTTT) and (|T5|) using different xc potentials. 

III. THREE-DIMENSIONAL TWO-BAND MODEL 

In general, it is very difficult to find the solution of the density matrix equation (TTTj) 
and one needs to make some approximations. For simplicity, we shall consider the optical 
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absorption spectra of systems composed of one valence and one conduction band: l\ = v, 
I2 = c. This approximation can be used if one assumes that the dominant optical transitions 
in the system take place from the highest occupied valence band to the lowest conduction 
band. Since the main aim of this paper is a proof of principle that TDDFT can describe 
excitonic effects in semiconductors triggered by short laser pulses, such an approximation 
is sufficient. However, for real systems one needs to take into account the band structure 
of the materials more accurately. Another simplification comes from the fact that in the 
dipole approximation for the external field, the optical transitions in the system take place 
with zero photon momentum. The coefficients c, defined in Eq. ( TTOl) . then depend on one 
momentum variable only: 

c^(t) = <wm 4 c q (t) = <wm (is) 

The problem is thus reduced to finding the density matrix of rank 2, with elements p™(t), 
p^{t), Pk'(^) an d Pk C (^)' which are functions of momentum and time. The elements p™(t) 
and p^it) describe the occupancy of the valence and conduction band states, and p^(t) and 
p^(t) describe the polarization in the system. These four elements are not independent. 
Particle conservation number requires p™{t) + p^{t) = 1, and by definition, p^(t) = p^*{t). 

From Eq. ffTT]) . one can obtain the following system of equations for two independent 
components p™(t) and p£ c (£): 

j t Pl v {t)= - 2Im[(E(t)dZ + VZl + V™ k )pi c (t)}, (19) 

§- t Pl c (t)= - •K-4]pr(*)-W(*)-pT(*M*)dk 

Wit) - pl v (t)}(v^(t) + v^(t)) - i\ygi(t) + v:ut) - vfii(t) - v x ^(t)]pi c (t). 

(20) 

This system represents the TDDFT version of the well-known SBEs,- which are used to 
study the optical properties of bulk semiconductors subject to external electric fields. Equa- 
tions ([191 . (120]) have a more general form, since here the Coulomb interaction effects are 
taken into account in principle exactly by means of the matrix elements of V xc [n](r,t). As 
mentioned in the Introduction, this is important in the case of small systems and sharp 
pulses, where the characteristic times in the system are shorter than the Coulomb scattering 
time, and where it is difficult to treat the Coulomb interaction effects properly within the 
BSE approach. 
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We shall study the solution of Eqs. (fT9l and (!20|) for the simple but instructive example 
of a two-band model on a cubic lattice. We assume that the solution of Eq. (j3J) with the 
Hamiltonian Eq. ([T]) gives the following simple dispersions: 

£ k = £ o + 2t v [cx)s(aok x ) + cos(a fc J/ ) + cos(a k z )}, (21) 

e^ = e C q — 2t c [cos(aok x ) + cos(aok y ) + cos(ao^)], (22) 

where ao is the lattice constant. For definiteness, we consider a model of solid hydrogen and 
assume that the difference of Eq and Eq can be set proportional to the first two hydrogen 
energy levels E = — l/[2asn 2 ],n = 1,2. In the last equation, a# = h 2 /me 2 = 1 is the 
Bohr radius, which is used as the length unit in this paper. It is assumed that the band 
parameters t v and t c are much smaller than Eq — i.e. the bands are non-overlapping. 
The valence and conduction band wave functions in the Bloch representation 

^(r)=e* k X(r), (23) 
^(r)=e* k X(r), (24) 

can be represented as linear combinations of Wannier functions w v (r) and w c (r). In fact, one 
can use for spatially periodic Bloch functions u^_(r) and u^(r) the following representation: 

ul( r ) = e lk(r - L V(r - L) (25) 

L 

and similar for «£.(!"), where L = ao(n xl n yi n z ) are the cell vectors. 

Unfortunately, it is difficult to find the exact expressions for the Wannier functions, 
and hence for the Bloch functions, even for such a simple three-dimensional system. For 
simplicity, we choose the Wannier functions to be equal to the Is and 2p° hydrogen wave 
functions: 

w v (r) = -^e- r , (26) 

\/7T 



w c (r) = ~^=e-^z. (27) 

This choice of the Wannier functions is an approximation, since the orthogonality condition 
on different sites J drw l *(r — L 1 )u? m (r — L 2 ) = ^ /m ^LiL 2 is violated. However, when the 
lattice parameter ao is much larger than the Bohr radius, as « ao, the requirement of 
orthogonality is satisfied with a very high precision, since the overlap between different sites 



becomes negligible. In our calculations, we shall use the lattice constant values do = 10 and 
20. 

In the next section, we shall study the solution of Eqs. (fl^|) and (1201) when an external 
short pulse field 

E(t) = E e- t2/r2 (28) 

with r ~ 1 — lOfs is applied. This will allow us to find the optical absorption spectrum A(u) 
of the system, which is defined by the real part of the ratio of the Fourier transforms of the 
total polarization 

P(u) = •(47r/ v ^)(eS - e%)\d™\ 2 al f J dte^tf® (29) 
and the pulse field E{oj): 

A{u) = -2Re[P(uj)/E(uj)]. (30) 



IV. RESULTS AND DISCUSSION 



A. Hartree-Fock 



Before proceeding with the solution of the TDDFT equations ( JT9l) . ( |20l) for different xc 
potentials, we briefly discuss the case of the Hartree-Fock potential for the two-band model 
presented above. Hartree-Fock is the most widely used approximation in the SBE theory^ 
for the study of optical absorption spectra in bulk systems and het erostructures subject to 
both smooth and pulsed external fields (see, for example, Refs. 



26 



23). 



In the time- dependent Hartree-Fock approximation, the total electron wave function sat- 
isfies the following equation: 



1 dt 



v 2 r 

-_-E(f)r + J dr 



(31) 



Applying the analysis presented in Section [TT] [see Eqs. (fTOfl - flTEj) ]. it is possible to show that 
Eq. (EH) is equivalent to the system of TDDFT density-matrix equations ( |T9l) . ( 1201) with 
corresponding matrix elements for the Hartree and Hartree-Fock potentials. Namely, the 
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Hartree potential matrix elements are: 

V&(*)= E 4 mnS (k,q)p q n W, (32) 

n,s=u,c;q 

where 

r r ib n *(r')'ib s (r') 

4 mns (k,q) =jdrj rf4'(r)C(r) | r _ r ?| • ( 33 ) 

From Eqs. (1321) . (133]) . one can find the following approximate expression for the Hartree 
matrix elements: 

V&(t) = 5 lm V(G - 0) + AlmnSa o / 7CT< n W + ( 34 ) 



n,s=v,c;ci 



where V(q) is the Fourier transform of the Coulomb potential l/|r — r'|, G is a reciprocal 
lattice vector and 

A imns = V(G) ! dvw l *{v)e- iGr w m {v) j dr'w n (v')e iGr 'w s (r') (35) 

G^O J J 



(details of a similar analysis can be found in Ref. l20f ) . In Eq. ( 134"]) . the dots correspond 
to the terms proportional to integrals over products of Wannier functions which reside on 
different sites [~ j drw l *(r)w m (r — L)]. These terms are negligible due to a small overlap of 
the Wannier functions on different sites [when as <C do, see Eqs. (1261) and ( 1271) ]. Moreover, 
the presence of oscillating functions exp[— zGr] and expfzGr'] under the integrals in the 
expression for the coefficients A ns , Eq. (|35|) . makes the terms proportional to A lmns much 
smaller than the first term in Eq. fl34|) . Therefore, the second term in Eq. fl34l) can be also 
neglected. Substitution of the remaining term 5 lm V(G 0) into Eqs. (fT9l) . (!20l) simply 
causes a constant shift of the overall potential and can be ignored. Therefore, the Hartree 
term does not contribute to the density-matrix equations ( [r9l) .( T20l) in this approximation, 
and we shall ignore it from now on. It is possible to show 20 that higher order terms in 
Eq. (1341) . which we neglect, lead to a small renormalization of the energy bands and the 
electric field time dependence. 

In a similar way, one can find an approximate expression for the matrix elements of the 
nonlocal Fock exchange potential: 

Vgh® = -al J j0j~ 3 V(k - q)pf (t). (36) 
11 



Substitution of the expression (136|) into Eqs. (|T9l) . (1201) leads to the following density-matrix 
equations in the Hartree-Fock approximation: 



21m 



+ a 3 



o 



(2k) 



^(k-q)(pr(t)-p^ c (t)M c W 



(37) 



- i[pm-pl V (t)) [E(t)<r k c -a, 



dq 



(38) 



These equations are equivalent to the standard SBEs.-) The solution of the system of 
Eqs. (1371) . (1381) in the presence of an external field E(t) allows one to calculate the opti- 
cal polarization in the system by means of Eq. ( J30l) . 

In Fig. [[J we present the frequency dependence of the optical absorption spectra at 
different values of the amplitude of the external field in Hartree-Fock approximation, when 
an external short-time electric pulse Eq. (1281) is applied. Here and in other Figures, we 




o.oo pv^'i 1 - 1 , i , 

-5.0 -2.5 0.0 2.5 5.0 

Frequency 



FIG. 1: Optical absorption spectra in Hartree-Fock approximation at different values of the valence 
bandwidth (divided by 7r 2 ). Here and in the Figures below, we consider external pulses Eq. (|28p of 
duration r. For this Figure, r = 0.05. The conduction bandwidth and the band gap are W c = 0.4-7T 2 
and uj g = 5. The amplitude of the electric field is Eq = 2.0. We consider the case of parabolic 
bands. 

introduce screening A = 1 in the Coulomb kernel l/(q 2 + A 2 ) and in order to reach a steady 
state more quickly we introduced a decoherence terms T(l + a\h\/ir)p lm on the right hand 
side of the density matrix equations. Here T = 0.2, a = 7.5. In this paper, in the numerical 
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results for the absorption spectra we set the prefactor (Stt / y/eb) (eq — EQ)\d cv \ 2 to be equal 1 
[see Eqs. ( |29|) and fl30|) ]. The absorption spectra demonstrate a pronounced excitonic feature 
with a shape and a peak position that depend on the external field amplitude parameters. 

The SBE approach was successfully applied to study optical properties of different bulk 
materials in the case of weak external fields, when linear response theory can be applied. 
Recently, this formalism was used to study the nonlinear response of different quantum-well 
systems,— >2£ and some other effects, like the sideband generation, Franz-Keldysh effect, four- 
wave mixing, and higher correlation effects, in particular, biexcitons etc (see, for example 
Ref. []). As discussed in the Introduction, the main shortcoming of the SBE approach for 
small (nano-) systems and short pulses is that the Coulomb interaction effects are treated 
the mean-field theory level. All the fluctuation effects are hidden in the inverse scattering 
(decoherence) time parameters T 1 ^ 1 , which are usually introduced into the system of the 
SBEs (1371) . (I38I) by means of the terms T^p 1 ^ 1 . However, the characteristic times in this 
case are shorter than the Coulomb scattering time; therefore, fluctuation effects are very 
important and can not be neglected. It is extremely difficult to include the higher Coulomb 
terms in the SBEs in a controllable way, since in this case one needs to consider additional 
equations for higher order correlation functions. The TDDFT approach allows one to treat 
the Coulomb effects in principle exactly, which should make this approach favorable in the 
case of small systems and short-time pulses, provided one can find a suitable xc functional 
to account for correlations. In the following two Subsections, we shall consider the optical 
properties of the two-band insulator by means of the TDDFT formalism developed in Section 
IITI We limit ourselves here to exchange-only functionals. 

B. Adiabatic LDA and GGA 

The LDA is a standard approximation used in DFT to study the ground state proper- 
ties of many materials. In the time-dependent case, there exists a generalization of this 
approximation - the adiabatic LDA (ALDA). In the adiabatic approximation it is assumed 
that the xc potential at time t depends on the time-dependent particle density n(r, t) at 
the same time t only. In other words, all memory effects are neglected, and the TDDFT 
equations remain local in time, which allows one to solve the problem numerically relatively 
easily. In ALDA, the x-only potential matrix elements have the following momentum and 
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time dependence: 

V£ DAk (t) = - (-) 1/3 I dv^n^t)^). (39) 

V 71 "/ J cell 

Our numerical analysis confirms the previous observation (see, e.g. Ref. Q) that the 
optical absorption spectra do not demonstrate excitonic features with the ALDA exchange 
potential (|39|) at various values of the model parameters (Fig. [2]). Actually, there is a very 
wide peak in the absorption spectra at rather large values of the quasiparticle masses (or 
very narrow energy bands); however, this wide absorption spectra weight cannot be related 
to a well-defined excitonic energy value, and hence it does not correspond to an exciton 
formation. For details of the numerical calculations in ALDA and other potentials used in 
this paper (see below), we refer the reader to Appendix lAl 




FIG. 2: Absorption spectra in ALDA at different values of the electric field amplitude and the free 
electron bands W v = 0.15, W c = 0.06. The other parameters are uo g = 0.375, clq = 20, t = 10, L = 
0.0025 and a = 7.5. 

Similar results can be obtained for standard adiabatic GGAs, for example the van 
Leeuwen-Baerends^ (LB) potential. In the LB approximation, the xc potential has the 
following form: 

V LB (r,t) = V^t) - finWfrt) f , x = ^ffi'?' (40) 

1 + 3/3x ln(x + Vx 2 + 1) nV 6 (r,t) 

where (3 is a parameter, which depends on the material. This potential has an important 
feature for finite systems, namely, a correct (~ 1/r) asymptotic behavior at long distances 
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(The ALDA potential decreases exponentially at r — > oo). In spite of this, no excitonic 
features are produced. 

Thus, we have seen that the standard adiabatic LDA+GGAs for the xc potentials does 
not allow one to reproduce the excitonic features in the optical absorption spectra within 
the TDDFT formalism developed above. One needs to go beyond this approximation and 
to find another class of simple TDDFT potentials which demonstrate excitonic peaks in 
the absorption spectra. As shown by Kim and Gorling^^ such potentials exist. The 
main requirement is a 1/q 2 singularity in the exchange energy kernel. In the following 
Subsection, we shall consider the optical absorption spectra for classes of functional where 
the underlying exchange energy has such a singularity. 



C. The Slater and KLI potentials 

As shown above, the ALDA cannot describe excitonic effects in the optical absorption 
spectra. A fundamental shortcoming of the LDA approximation is that it contains a self- 
interaction error. In order to reduce this error, several approximations for self-interaction 
corrected (SIC) potentials, were proposed. Probably the most often used SIC potential is 
due to Perdew and Zunger.— However, in their scheme one must deal with orbital-dependent 
potentials, which makes the calculations difficult, especially in the time-dependent case. To 
make a computational procedure simpler, the method of the optimized effective potential 
(OEP) was proposed (for overview and references, see for example Ref.— ). In this approx- 
imation, all the wave functions for different orbitals satisfy a Schroedinger equation with 
a common, orbital-independent, xc potenti al v^ p (r). A time-dependent generalization of 
the OEP approximation was given in Refs. 
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31 



Unfortunately, this method is computa- 
tionally very demanding in the time- dependent case. Krieger, Li, and Iafrate^ proposed 
a simplified method (the KLI scheme) for the OEP in the equilibrium exact-exchange case 
with v°^ p (y) depending explicitly on the orbital functions ipj a . In the equilibrium case, the 
KLI potential is defined by the following integral equation: 



E^y{ M ^W + / ^IMOI 2 [«^(rO-^r')]}, (41) 



were the orbital-dependent potentials u XC j a (r) are obtained from the xc energy E xc as follows: 



u 



XCJ 



ff (r) = [f 3a ^ a (v)}-'5E xc [{ Vja }}/5 Vja (v)- (42) 
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This approximation was generalized to the time dependent case in the adiabatic approxi- 
mation (n(r) — ► n(r,t)) in Ref. |31. 

In the linear regime, the optical absorption spectra of insulators show excitonic peaks, 
when an OEP approximation with the exchange energy with a kernel that contains a 1/g 2 
singularity is used.— ^ However, it is much more difficult to use this approach for short and 
strong pulses, when one needs to go beyond linear response. 

Here, we consider simplified OEPs with a Coulomb singularity in the exchange energy 
kernel, which are approximate solutions of Eqs. (14T]) . flUD , and solve the problem in the time- 
domain by using the TDDFT formalism developed in the previous Sections. This approach 
can be applied for external fields of arbitrary strength and duration. 

We consider the KLI potential and the Slater potential^, which can be obtained if one 
neglects the orbital-dependent constants on the right hand side of Eq. (1411) : 

n ja (r,t) 



Slater/ 



n a (r,t) 



■u 



xcja(r,t) • 



(43) 



We use the Fock exchange energy for our two-band model in Eq. (1421) . Optical absorption 
spectra obtained with the Slater and KLI exchange potentials at different values of the model 
parameters are presented in Figs. [3] and HI 
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FIG. 3: Optical absorption spectra in Slater approximation at different values of the valence 
bandwidth (divided by 7r 2 ). The other parameters are W c = 0.0257T 2 , uj g = 0.2, a = 10, r = 4, E = 
0.0166, T = 0.012 and a = 7.5. 

As can be seen from these Figures, the absorption spectra demonstrate a significant 
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excitonic peak. In particular, the excitonic peak amplitude and the exciton binding energy 
decrease with increasing valence bandwidth in the Slater case (Fig. [3]). This effect takes 
place because the increase of the bandwidth is equivalent to a reduction of the effective mass 
of electron-hole pairs. A similar behavior can be found in the case of the KLI potential. 
Fig. H] shows the KLI absorption spectrum at different values of the electric field. Here, the 
spectrum demonstrates nonlinear effects when the pulse amplitude is large. 

There is a significant quantitative difference between different xc potentials. In partic- 
ular, we find that the excitonic features become more pronounced as one passes from KLI 
to Slater. The excitonic effects are defined mainly by the ratio of the Coulomb interaction 
energy to the typical energy of the free system. Since in our case the valence and conduc- 
tion bands and the gap have the same order of magnitude, we can consider the ratio of the 
Coulomb energy to the valence bandwidth. In the Slater approximation, the amplitude of 
the Coulomb interaction energy appears to be bigger compared to KLI. In fact, as follows 
from the Appendix, the Coulomb interaction matrix elements are much larger in the Slater 
approximation (see Eqs. (lAlOf) . (1A12I) and ( 1A13I) ). In the strongly localized case, consid- 
ered in this paper, the Slater results are closer to HF (Fig. [5]). For less localized electrons, 
when ao ~ 1 and the momentum dependence of the xc potential matrix elements can not 
be neglected, we expect to get closer agreement between Slater and KLI. Another impor- 
tant conclusion from our numerical results is a decrease of the excitonic binding energy in 
Slater and KLI compared to HF. In fact, it is known that the HF app roximation gives an 



overestimated binding energy of excitons (for a discussion, see Ref. l33l ). 

Thus, we have shown that TDDFT can describe excitonic effects in the optical absorption 
spectra triggered by short-time laser pulses, when the exchange energy kernel contains a 
long-range Coulomb singularity. 



V. CONCLUSIONS 



In this paper, a density-matrix TDDFT formalism to describe ultrafast processes in 
semiconductor structures has been developed. We have shown that the corresponding system 
of equations for the matrix elements is a generalization of the SBEs, since it allows one to 
take into account Coulomb effects in principle exactly. This is very important in the case 
of short pulses and when the system is very small, such as for nanostructures. Another 
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Frequency 

FIG. 4: Optical absorption spectra in KLI approximation at different values of the electric field. 
The other parameters are W v = 0.0416vr 2 , VF C = 0.0208tt 2 , co g = 0.166, a = 10, r = 4.8, F = 0.012, 
and q = 7.5. 

important advantage of this formalism is that it can be applied directly in the time domain, 
which makes the calculation much faster. As an example, we studied the optical absorption 
spectra of a two-band model bulk insulator. We have shown that the optical absorption 
spectra include excitonic effects when the xc energy kernel contains a 1/q 2 singularity. The 
exciton binding energy within the Slater and KLI approaches is much smaller compared to 
HF, which often gives too large value for the binding energy. 

This TDDFT approach should be useful for studying ultrafast processes in real semicon- 
ductor and polymer nanostructures, for example for short-time laser pulse experiments. In 
particular, it gives access to a variety of nonlinear effects, which will be the subject of future 
studies. 
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FIG. 5: Optical absorption spectra in different approximations at the same values of the model 
parameters. The valence bandwidth is equal to 0.05-7T 2 and the other parameters are given in the 
caption to Fig. [3j The Hartree-Fock result is obtained using a Bloch function approximation for 
the xc potential matrix elements (see Appendix). 

APPENDIX A: TECHNICAL DETAILS 

In this Appendix, we present some useful approximations and technical details of the 
numerical solution of the TDDFT equations (fl"9l . (!20|) . Since, generally speaking, evaluation 
of the xc matrix elements beyond LDA requires six-dimensional spatial integration, it is 
necessary to make some simplifications to speed up the calculations. We shall use the 
Wannier and Bloch representation ( |23l) -( l27l) for the wave functions. In this case, the electron 
density can be expressed in the following form: 

n(r,t) = 2 P l ™(t)e- ih(l w m *(r - L)w l (v). (Al) 

Z,m,L,q 

Since the wave functions w m {r) and w l (r) are strongly localized, the main contribution to 
expression (lAip is given by the term L = 0. Therefore, the electron density is approximately 

n(r,t)^2^p^(tV m *(rV(r), (A2) 

l,m 

where 

P^)=5>?^)- (A3) 

q 
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It turns out that Eq. (1A2j) can be obtained by another approximation for the wave 
functions, which follows from Eqs. fl23]) - fl25]) . Namely, 

€(^e'X(r), (A4) 
^W^XW. (A5) 

We choose the wave vector ko to be zero. This approximation is good when one considers 
excitation processes around the direct gap. In this case, 

<(r)^e* k W(r), (A6) 
V£(r) ~ e ikr w c (r). (A7) 

By using this approximation one gets the following expression for the time-dependent matrix 
elements of the LDA and Slater potentials: 



V&*k(*) * - ^y /3 J^dvw l *(T)w m (r) 



1/3 

w a *(r)w\r) ) , (A8) 



V l s TateAt) ~ - Y. E V <P - *)A lmnsn *(t)p™(t)p™(t), (A9) 

n,s,n,s p,q 

where 

A imnsns (t) „ f ^ w 1 * (r)w m (r)w n * (v)w s (r)w n * (v)w~ s (r) 

U J cea n{v,t) 1 ] 

and the density n(r, t) is given in Eq. (1A2|) . 

Finally, by using the Wannier wave functions (125|) - (|27|) and making an approximation 
similar to that used in the derivation of Eqs. (IA8I) - (IA9I) . one gets the following equation for 
the KLI matrix elements, which follows from Eq. f j4Tl ): 



n,s,n,s,p n,fi,s,s,p,q 

(All) 

where Vsiater-pif) is defined in Eqs. (IA9[) and flAlOl) . and 

B l r s [t) = f ^(T)w"(T)wy)w*(v) ^ (A12) 



cell 



n(r, t) 
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Bl ~ {t) ^_f rfr ^(r)^(r)^(r)^(r)^(r)^(r)^ wa 
J ceii n[r,t) ^ 

, u;'*(r)w m (r)w n *(r)w s (r)w"'*(r)w 5 (r)M; , '*(r)w , '(r) 

i^Frt • ( 3) 

cell n \ r j l ) 

It is convenient to solve Eq. flAllj) by iteration. Therefore, we have reduced the six- 
dimensional space integral to a three-dimensional one. For numerical integration, we divide 
the space interval into 120 parts for every direction. In the case of momentum integration, 
we divide the interval into 100 parts. As it was mentioned in Section II V Al we use a deco- 
herence factor T in equations in order to reach a steady state faster. We have found that 
the position of the excitonic peak does not depend on value of T, but its height decreases 
with T increasing. Also, it was found that in the Slater and KLI cases the total absorption 
conservation law is satisfied with a higher precision when T increases. Due to finite values of 
the time interval and time step, absorption A(u) acquires unphysical finite weights at large 
values of \uj — uj 9 \. These weights disappear when Y increases. 

For completeness, let us present approximate matrix elements for the Hartree-Fock po- 
tential derived within the same approximation. These expressions are necessary to make a 
comparison of the absorption spectra using different potentials (Fig. [5]) . The Hartree-Fock 
matrix elements have the following structure: 

VwicCO ^ - E E V ^ k - P)A lmns C(t), (A14) 

n,s p,q 

where 



A lmns ~ / drw l *{r)w m {r)w n *(r)w s (r). (A15) 

J cell 

Following the discussion in Sec. IIV Al the quantities ( 1A15I) can be approximated as A lmns ~ 
5 ls 5 mn . However, in approximation (1A4I) . (1A5I) these elements become 

A lmns _ 5 mnjm f drw l * ( T )w 9 (r) ~ ^Vi™, (A16) 
J cell 

where I m ,m = v,c, is the main contribution of the function |w m (r)| 2 to the integral 
Eq. (TAT51) . Its value can be chosen to be equal to the maximum value of |w m (r)| 2 , i.e. 
1/tt for the valence band function and l/(167re 2 ) for the conduction band function. In the 
last case, there is an additional factor 1/2, which comes from the azimuthal angle inte- 
gration. The factor 7r/6 must be added in Eq. (IA15I) . since in Slater and KLI there is an 
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additional momentum integration over the Brillouin zone, and the integration is performed 
over a sphere of the radius 7r instead of the cubic Brillouin zone. 

The most important property of expressions (1A8[) - (1A9j) and (lAlip is the fact that they are 
momentum-independent. Therefore, as follows from Eqs. (fT§|) and f[2"Uj) . the matrix elements 
Pq (t) are functions of the energy 



Ek = - [cos k x + cos k y + cos k z 



(A17) 



In order to reduce the dimensionality of the integrals it is convenient to introduce the 
corresponding three-dimensional density of states 

d 3 p 



Die] 



(2tt) 3 



8[e — (1/3) (cos q x + cosq y + cosq z )} 
71 9[1 — (3e — cosp x — cospy) 2 } 



3 r , r , en - 

— dp x dp y —= 

7T 3 Jo Jo \/l 



, 7 . (A18) 

a/1 — (3e — cosp x — cospy) 2 

where the energy e changes from —1 to 1. The energy dependence of the density of states 

is shown in Fig. 




-0.5 0.0 0.5 

energy 



FIG. 6: The three-dimensional density of states Eq. (IA18j) as a function of energy. 



Since the xc matrix elements depend on integrals which have the structure 



d 3 p f d 3 q 



;V(p-q)F(p^ jP ^ 



(2tt) 3 J (2tt) 3 

it is also convenient to introduce a two-energy density of states 

d 3 p f d 3 q 



(A19) 



D 2 (e,e) 



(2vr) 3 J (2tt 



5[e — (1/3) (cos ^ + cospy + cosp^)] 
xS[e — (l/3)(cosg x + cosq y + cosq z )]V(p — q). (A20) 
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Performing the integration over p z and q z in Eq. (1A20I) . one finds 




VI + 5 - x 2 y/1 + 5 - y 2 



6[l-x 2 } 6[l-y 2 } 



1 



x 



(pil — qy) 2 + (cos 1 x — cos 1 y) 2 + \/X 



2 ' 



(A21) 



where x = 3e — cosp x — cosp y ,y = 3e — cosq x — cosq y and py = (p x ,P y ,0). In Eq. (1A21I) . 
we have introduced "screening" parameters 5 and A, since in the three-dimensional case 
the function D2(e,e) has a logarithmic singularity at e — e. In fact, it is possible to show 
analytically that in the case of a square dispersion D 2 (e — > e) ~ \n(l/\,/e — \fi\). We 
apply this approximation in the LDA case (Fig. (T5])). In other Figures, we use the parabolic 
band approximation for sake of simplicity. This approximation gives results close to those 
obtained by using the densities of states Eqs. (1A18|) and (1A21j) . when one uses a decoherence 
factor with a momentum cutoff (see Section TlV Al) . In fact, in this case low momenta give a 
dominant contribution to the absorption processes. 
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